clear
%constant
global gamma Nx eps length dx dt
gamma = 1.4;
length = 2;
time = 0.14;
eps = 1e-6;
% discretization
Nx = 261;
Nt = 350;
dx = length/Nx;
dt = time/Nt;
rho_array = zeros(1,Nx);
rho_u_array = zeros(1,Nx);
E_array = zeros(1,Nx);
%initialization
for i = 1:1:Nx
    if i <= floor(Nx/2)
        rho_array(i) = 0.445;
        rho_u_array(i) = 0.311;
        E_array(i) = 8.928;
    else
        rho_array(i) = 0.5;
        rho_u_array(i) = 0;
        E_array(i) = 1.4274;
    end
end
v_array = [rho_array; rho_u_array; E_array];

%time stepping
for i = 1:1:Nt
    v_array = v_array + time_step(v_array);
    i
end

rho_array = v_array(1,:);
rho_u_array = v_array(2,:);
E_array = v_array(3,:);
plot(rho_array,'-o','markersize',3);
hold on;
plot(rho_u_array,'-o','markersize',3);
hold on;
plot(E_array,'-o','markersize',3);

function dv_array = time_step(v_array)
    global Nx gamma eps length dx dt
    rho_array = v_array(1,:);
    rho_u_array = v_array(2,:);
    E_array = v_array(3,:);
    %auxiliary parameters p,a, 
    u_array = rho_u_array./rho_array;
    p_array = (gamma-1)*(E_array-1/2*rho_array.*u_array.^2);
    a_array = speed(rho_array,p_array);
    lambda_array = [u_array; u_array-a_array; u_array+a_array];

    lambda_pos = (lambda_array + sqrt(lambda_array.^2 + eps))/2;
    lambda_neg = (lambda_array - sqrt(lambda_array.^2 + eps))/2;

    fpos_array = SW_split(lambda_pos, u_array, rho_array, p_array, a_array);
    fneg_array = SW_split(lambda_neg, u_array, rho_array, p_array, a_array);
    flux_lim = limited_flux(u_array,rho_array,E_array,a_array);
    for i = 1:1:Nx-1
        flux(:,i) = fpos_array(:,i) + fneg_array(:,i+1) + flux_lim(:,i);
    end
    % Upwind scheme
    
    for i = 2:1:Nx-1
        drho_array(i) = - (dt/dx)*(flux(1,i) - flux(1,i-1));
        drho_u_array(i) = - (dt/dx)*(flux(2,i) - flux(2,i-1));
        dE_array(i) = - (dt/dx)*(flux(3,i) - flux(3,i-1));
    end
    % boundary points
    drho_array(1) = drho_array(2);
    drho_u_array(1) = drho_u_array(2);
    dE_array(1) = dE_array(2);
    drho_array(Nx) = drho_array(Nx-1);
    drho_u_array(Nx) =drho_u_array(Nx-1);
    dE_array(Nx) = dE_array(Nx-1);
    dv_array = [drho_array; drho_u_array; dE_array];
end

%speed of sound
function a = speed(rho, p)
    global gamma Nx 
    for i = 1:1:Nx
        a(i) = sqrt(gamma*p(i)/rho(i));
    end
end
%linear transform of u
function f = SW_split(lambda, u, rho, p, a)
    global gamma Nx 
    for i = 1:1:Nx
        w(i) = (3-gamma)*(lambda(2,i)+lambda(3,i))*a(i)^2/(2*(gamma-1));
        f(1,i) = 2*(gamma-1)*lambda(1,i) + lambda(2,i) + lambda(3,i);
        f(2,i) = 2*(gamma-1)*lambda(1,i)*u(i) + lambda(2,i)*(u(i)-a(i)) +lambda(3,i)*(u(i)+a(i));
        f(3,i) = (gamma-1)*lambda(1,i)*u(i)^2 + lambda(2,i)/2*(u(i)-a(i))^2 + lambda(3,i)/2*(u(i)+a(i))^2 + w(i);
        f(:,i) = f(:,i)*(rho(i)/(2*gamma));
    end
end

function fl = limited_flux(u, rho, E, a)
    global gamma Nx dx dt
    fl = zeros(3,Nx-1);
    alpha_array = zeros(3,Nx-1);
    lambda_array = zeros(3,Nx-1);
    Right_array = zeros(3,3,Nx-1);
    for i = 1:1:Nx-1
        u_hat = (u(i)*rho(i)^0.5+u(i+1)*rho(i+1)^0.5)/(rho(i)^0.5+rho(i+1)^0.5);
        H_hat = ((a(i)^2/(gamma-1)+0.5*u(i)^2)*rho(i)^0.5+(a(i+1)^2/(gamma-1)+0.5*u(i+1)^2)*rho(i+1)^0.5)/(rho(i)^0.5+rho(i+1)^0.5);
        a_hat = ((gamma-1)*(H_hat-0.5*u_hat^2))^0.5;
        Right_mat = [1, 1, 1;
                     u_hat-a_hat, u_hat, u_hat+a_hat;
                     H_hat-u_hat*a_hat, 0.5*u_hat^2, H_hat+u_hat*a_hat];   
    
        Delta_U = [rho(i+1)-rho(i); rho(i+1)*u(i+1)-rho(i)*u(i); E(i+1)-E(i)];
        alpha = Right_mat \ Delta_U;
        lambda = [u_hat-a_hat; u_hat; u_hat+a_hat];
        alpha_array(:,i) = alpha;
        lambda_array(:,i) = lambda;  
        Right_array(:,:,i) = Right_mat;
    end
    for i = 2:1:Nx-2
        lambda = lambda_array(:,i);
        for p = 1:1:3
            if lambda(p) > 0 && abs(alpha_array(p,i))>1e-6
                theta = alpha_array(p,i-1)/alpha_array(p,i);
            elseif lambda(p) < 0 && abs(alpha_array(p,i))>1e-6
                theta = alpha_array(p,i+1)/alpha_array(p,i);
            else
                theta = 1;
            end
            alpha_l = alpha_array(p,i)*limiter_function(theta);
            fl(:,i) = fl(:,i) + 0.5*abs(lambda(p))*(1-dt/dx*abs(lambda(p)))*alpha_l*Right_array(:,p,i);
        end
    end
end

% % Van leer limiter
% function phi = limiter_function(theta)
%     phi = (theta + abs(theta))/(1+abs(theta));
% end

% % superbee limiter
% function phi = limiter_function(theta)
%     phi = max(max(0, min(1,2*theta)),min(2,theta));
% end

% MC limiter
function phi = limiter_function(theta)
    phi = max(0, min(min(2,2*theta),(1+theta)/2));
end

% % minmod limiter
% function phi = limiter_function(theta)
%    if abs(theta)>0
%        if theta > 1
%            phi = 1;
%        else
%            phi = theta;
%        end
%    else
%        phi = 0;
%    end
% end
     